include("program.jl")
using Plotly
importall DoubleIntegralUtils
importall DoubleIntegralNewtonCotes
importall DoubleIntegralGauss
importall DoubleIntegralPadua
function draw_functions(a,b,c,d,precision, function_list)
function get_scatter(f,a,b,c,d,precision)
domain3d = [(x,y) for x in linspace(a,b,precision) for y in linspace(c,d,precision)]
return Plotly.scatter3d(x=[t[1] for t in domain3d],y=[t[2] for t in domain3d],z=[(f(t[1], t[2])) for t in domain3d], mode="lines")
end
Plotly.plot([get_scatter(func,a,b,c,d,precision) for func in function_list])
end
function test_function(f,a,b,c,d,exact_value,n)
ns = [2:1:n;]
Plotly.plot([Plotly.scatter(x=ns, y=[abs(composite_double_trapezoid(f,a,b,c,d,i,i)-exact_value)/abs(exact_value) for i in ns], name="Composite Trapezoid Method Error"),
Plotly.scatter(x=ns[1:2:end], y=[abs(composite_double_simpson(f,a,b,c,d,i,i)-exact_value)/abs(exact_value) for i in ns[1:2:end] ], name="Composite Simpson Method Error"),
Plotly.scatter(x=ns[1:1:min(n-1,14)], y=[abs(double_n_by_m_order(i,i,f,a,b,c,d)-exact_value)/abs(exact_value) for i in ns[1:1:min(n-1,14)]], name="Newton-Cotes Cubature Error"),
Plotly.scatter(x=ns, y=[abs(double_gauss_chebyshev((x,y) -> (move_function(f,a,b,c,d)(x,y))*sqrt(1-x^2)*sqrt(1-y^2),i,i)-exact_value)/abs(exact_value) for i in ns], name="Gauss Cubature Error"),
Plotly.scatter(x=ns, y=[abs(padua(move_function(f,a,b,c,d), i)-exact_value)/abs(exact_value) for i in ns], name="Padua Cubature Error")])
end
function test_exact_digits(f,a,b,c,d,exact_value,n)
ns = [2:1:n;]
Plotly.plot([Plotly.scatter(x=ns, y=[exact_digits(composite_double_trapezoid(f,a,b,c,d,i,i),exact_value) for i in ns], name="Composite Trapezoid Method Exact Digits"),
Plotly.scatter(x=ns[1:2:end], y=[exact_digits(composite_double_simpson(f,a,b,c,d,i,i),exact_value) for i in ns[1:2:end] ], name="Composite Simpson Method Exact Digits"),
Plotly.scatter(x=ns[1:1:min(n-1,14)], y=[exact_digits(double_n_by_m_order(i,i,f,a,b,c,d),exact_value) for i in ns[1:1:min(n-1,14)]], name="Newton-Cotes Cubature Exact Digits"),
Plotly.scatter(x=ns, y=[exact_digits(double_gauss_chebyshev((x,y) -> (move_function(f,a,b,c,d)(x,y))*sqrt(1-x^2)*sqrt(1-y^2),i,i),exact_value) for i in ns], name="Gauss Cubature Exact Digits"),
Plotly.scatter(x=ns, y=[exact_digits(padua(move_function(f,a,b,c,d), i),exact_value) for i in ns], name="Padua Cubature Exact Digits")])
end
Ze względu na limity dotyczące rozmiaru plików, zachęcam do wykorzystania zakomentowanego kodu.
Funkcja move_function(f,a,b,c,d) zwraca taką funkcję g, że $ \int_a^b \int_c^d f(x,y) dydx = \int_{-1}^{1} \int_{-1}^{1} g(x,y) dydx $.
f(x,y) = 1/(1+x^2+y^2)
g = move_function(f,0,1,0,1)
gf = (x,y) -> g(x,y)*sqrt(1-x^2)*sqrt(1-y^2)
@show composite_double_trapezoid(f,0,1,0,1,10,10)
@show composite_double_simpson(f,0,1,0,1,10,10)
@show Float64(double_n_by_m_order(10,10,f,0,1,0,1))
@show double_gauss_chebyshev(gf, 10,10)
@show padua(g, 10)
f(x,y) = 1/(x^2 + y^2 + 1)
draw_functions(0,1,0,1,100, [f])
exact_int = 0.639510351870311001962693085427323679
test_function(f,0,1,0,1,exact_int, 30)
test_exact_digits(f,0,1,0,1,exact_int, 30)
f1(x,y) = 0.75 * exp(-(9*x-2)^2/4 - (9*y-2)^2/4)
f2(x,y) = 0.75 * exp(-(9*x+1)^2/49 - (9*y+1)/10)
f3(x,y) = 0.50 * exp(-(9*x-7)^2/4 - (9*y-3)^2/4)
f4(x,y) = -0.2 * exp(-(9*x-4)^2 - (9*y-7)^2)
f(x,y) = f1(x,y) + f2(x,y) + f3(x,y) + f4(x,y)
exact_int = 0.40696958949155611
draw_functions(-1,1,-1,1,100, [move_function(f,0,1,0,1)])
test_function(f,0,1,0,1, exact_int, 30)
test_exact_digits(f,0,1,0,1, exact_int, 30)
Możemy sprawdzić, że $\int_{-1}^1 \int_{-1}^1 \sin(x+y+1) dydx = 4sin^3(1)$.
f(x,y) = sin(x+y+1)
exact_int = 4*sin(1)^3
#draw_functions(-1,1,-1,1,100, [f])
#test_function(f,-1,1,-1,1,exact_int, 30)
#test_exact_digits(f,-1,1,-1,1,exact_int, 30)
f(x,y) = exp(x+y+1)
exact_int = e*(e-1/e)^2
#draw_functions(-1,1,-1,1,100, [f])
#test_function(f,-1,1,-1,1,exact_int, 30)
#test_exact_digits(f,-1,1,-1,1,exact_int, 30)
f(x,y) = sqrt((x+1)*(y+1))
exact_int = 32/9
draw_functions(-1,1,-1,1,100, [f])
test_function(f,-1,1,-1,1,exact_int, 12)
test_exact_digits(f,-1,1,-1,1,exact_int, 30)
f(x,y) = 1/(1 + 25x^2 + 25y^2)
exact_int = 0.4361934108122879
test_function(f, -1, 1, -1, 1, exact_int, 15)